/*
 * Multiphase.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "Multiphase.h"

// C++ includes
#include <algorithm>
#include <iomanip>
#include <tuple>

Multiphase::Multiphase()
{}

Multiphase::Multiphase(const vector<Phase>& phases) :
phases(phases)
{
	// Loop over all the phases to set the sequence of species 
	BOOST_FOREACH(const Phase& phase, phases)
	{
		// The species in the current phase
		auto speciesPhase = phase.GetSpecies();
		
		species.insert(species.end(), speciesPhase.begin(), speciesPhase.end());
	}
	
	speciesMap = SpeciesMap(species);
}

const vector<Phase>& Multiphase::GetPhases() const
{
	return phases;
}

const vector<string>& Multiphase::GetSpecies() const
{
	return species;
}

const SpeciesMap& Multiphase::GetSpeciesMap() const
{
	return speciesMap;
}

template<class PhaseFunction> 
const VectorXd GenericPhaseComputation(double T, double P, const VectorXd& n, const vector<Phase>& phases, const PhaseFunction& phaseFunction)
{
	// The activity-like data of the species to be computed
	VectorXd a(n.rows());
	
	// The current number of processed species, to be update at the end of every loop (see below)
	uint nprocessed = 0;
	
	// Loop over all the phases
	BOOST_FOREACH(const Phase& phase, phases)
	{
		// The number of species in the current phase
		const uint nspecies = phase.GetSpecies().size();
		
		// The vector of number of moles of the species of the current phase
		const VectorXd np = n.segment(nprocessed, nspecies);
		
		// Compute the activity-like data of the species in the current phase using the template parameter phaseFunction
		const VectorXd ap = phaseFunction(phase, T, P, np);
		
		// Transfer the activity-like data of the current phase to the vector of activity-like data a (beginning at position given by nprocessed)
		a.segment(nprocessed, nspecies) = ap;
		
		// Update the number of processed species
		nprocessed += nspecies;
	}

	return a;
}

const VectorXd Multiphase::Concentrations(const VectorXd& n) const
{
	auto concentrationsFunction = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.Concentrations(n); };
	
	return GenericPhaseComputation(0, 0, n, phases, concentrationsFunction); // T = 0, P = 0 because their values are meaningless for the calculation of the concentrations
}
	
const VectorXd Multiphase::ActivityCoefficients(double T, double P, const VectorXd& n) const
{
	auto activityCoefficientsFunction = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.ActivityCoefficients(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activityCoefficientsFunction);
}
	
const VectorXd Multiphase::Activities(double T, double P, const VectorXd& n) const
{
	auto activitiesFunction = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.Activities(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activitiesFunction);
}

void Multiphase::PrintState(double T, double P, const VectorXd& n) const
{
	VectorXd c = Concentrations(n);
	VectorXd g = ActivityCoefficients(T, P, n);
	VectorXd a = Activities(T, P, n);
	
	typedef std::tuple<string,double,double,double,double> Entry;
	
	vector<Entry> data(species.size());
	
	for(unsigned i = 0; i < species.size(); ++i)
		data[i] = Entry{species[i], n[i], a[i], g[i], c[i]};
	
	std::sort(data.begin(), data.end(), [](const Entry& e1, const Entry& e2) { return std::get<1>(e1) > std::get<1>(e2); });
	
	cout << setprecision(6);
	
	cout << setw(25) << left << "Species";
	cout << setw(25) << left << "Moles";
	cout << setw(25) << left << "Activities";
	cout << setw(25) << left << "ActivityCoefficients";
	cout << setw(25) << left << "Concentrations" << endl << endl;
	
	BOOST_FOREACH(const Entry& e, data)
	{
		cout << setw(25) << left << std::get<0>(e);
		cout << setw(25) << left << std::get<1>(e);
		cout << setw(25) << left << std::get<2>(e);
		cout << setw(25) << left << std::get<3>(e);
		cout << setw(25) << left << std::get<4>(e) << endl;
	}
}

ostream& operator<<(ostream& out, const Multiphase& multiphase)
{
	const SpeciesMap& speciesMap = multiphase.GetSpeciesMap();
	
	BOOST_FOREACH(const Phase& phase, multiphase.GetPhases())
	{
		out << "Phase: " << phase.GetName() << endl << '\t';
		
		BOOST_FOREACH(const string& s, phase.GetSpecies())
			out << "(" << speciesMap[s] << ") " << s << ", ";
		
		out << endl << endl;
	}
	
	out << "Species: " << endl << '\t';
	
	BOOST_FOREACH(const string& s, multiphase.GetSpecies())
		out << "(" << speciesMap[s] << ") " << s << ", ";
	
	return out;
}
